建筑物是城市空间重要的地理信息要素，从遥感数据中自动提取建筑物对于城市三维建模、灾害评估、数字地图更新等应用具有重要的作用[1-2]。建筑物的自动提取一直是摄影测量与遥感领域一个重要的研究方向，但是由于城市环境复杂，从遥感数据中自动提取建筑物仍然面临许多困难[3]。

近年来，国内外学者对建筑物提取进行了广泛的研究。根据所用数据的不同，可以分为基于影像数据[4-6]、基于LiDAR数据[7-10]，以及联合影像与LiDAR数据[11-15]3类方法。第1类方法主要利用影像的光谱信息。由于影像广泛存在的“同物异谱”、“同谱异物”现象，以及阴影、遮挡等现象，单纯利用影像提取的建筑物有较多的错误，难以实用化[16]。利用LiDAR数据提取建筑物主要利用点云的强度、回波以及几何信息，虽然自动化程度较高，但相较于影像数据，LiDAR数据的空间分辨率较低，提取的建筑物轮廓不够精确。

融合LiDAR点云与影像数据提取建筑物，可以充分利用LiDAR数据精确的三维信息和影像的光谱信息以及分辨率高的优势，二者的融合更加有利于自动提取高精度的建筑物信息。文献[11]提出一种面向对象的航空影像与LiDAR数据融合分类方法，但是该方法依赖于高精度的航空影像分割。文献[12]首先利用光谱信息分离出植被，然后利用高程信息从地面道路和建筑物中分离出建筑物。该方法单独利用光谱信息分离植被，利用高程信息分离建筑物，没有将两种数据进行有效融合。文献[13]采用支持向量机融合光谱与几何特征进行建筑物和植被提取，但是监督分类方法需要大量的标记数据。以上方法主要利用了影像的光谱信息来提高建筑物提取精度，而忽略了影像的高分辨率特征来优化建筑物边缘信息。文献[14]利用超高分辨率影像在LiDAR辅助下提取建筑物轮廓；文献[15]融合LiDAR与正射影像进行建筑物提取和边缘规则化。两种方法都是基于影像上的线特征提取进行边缘信息规则化，但是线特征提取通常不完整，易受阴影区域的影响。针对现有方法存在的问题，基于LiDAR点云与影像数据的特点，本文引入图割算法[17]，提出一种融合几何信息与颜色信息的建筑物提取方法，并利用前后景分割的思想来优化建筑物边缘。在本文方法中，一方面通过图割算法来融合几何与颜色信息，考虑邻域关系以保证提取结果邻域内的一致性；另一方面，通过前后景分割的方法基于影像的颜色信息来优化建筑物边缘，以利用影像分辨率高的优势。本文方法主要流程如图 1所示。

仪器误差或者航空激光点云获取过程中存在的飞鸟，使得获取的点云中存在粗差点，影响后续的DSM内插以及几何特征计算，首先利用基于统计分析的方法剔除粗差点[18]。在此基础上，采用文献[19]的方法将点云分为地面点和非地面点。然后利用全部点云内插DSM，利用地面点内插DTM，再利用DSM减去DTM获取nDSM。由于本文采用的图割优化过程在DSM格网上实现，因此提出一种保边缘特征的DSM内插方法，主要包括以下几步：

(1) 由点云平均密度估计DSM格网分辨率R。

(2) 含有LiDAR点的格网高程采用其包含的LiDAR点的最大高程值。

(3) 缺少LiDAR点的格网高程通过邻近格网计算。首先将当前格网周围8邻域内的有效格网(含有LiDAR点的格网)的高程按照1.5 m间隔统计直方图。当某一个区间内含有的格网数大于有效格网数的一半时，直接将这一区间格网高程的平均值赋予当前格网，否则利用8邻域格网中所有的LiDAR点，按照与格网中心距离反距离加权计算当前格网高程。

(4) 如果8邻域内没有有效格网，则继续扩大邻域范围，直到内插出有效高程。

为节省后续处理时间，将非地面点云投影到DSM格网上，生成一张非地面掩膜Mng以指示非地面区域，也就是建筑物提取候选区域。在Mng中，“0”指示地面区域，“1”指示非地面区域。

NDVI常用于区分植被和非植被区域，本文也将NDVI作为区分植被和建筑物的一个特征。但是阴影区域较暗的颜色会对NDVI的计算产生较大影响[20]，对此，首先采用文献[21]中的方法检测阴影区域，然后对阴影区域进行增强。具体做法是利用连通分析将阴影区域聚类，然后对每个聚类的对象分别进行阴影增强，即将影像的RGB转换到HSI颜色空间，然后按照式(1)对S和I通道进行拉伸，再转换到RGB空间，以实现对阴影区域的增强。

本文主要采用了基于点云的3个几何特征：点云平整度、法向量分布、基于nDSM影像的纹理一致性，以及基于正射影像的颜色特征NDVI。下面对采用的特征进行介绍：

通常情况下，建筑物由多个平面组成，而植被区域表面不规则，因此建筑物上的LiDAR点云具有局部平整性。本文采用协方差矩阵分析来计算点云这一特征[22]。

令NP=pi|i=1, 2, …, n表示非地面点，pi=(xi, yi, zi)为其中一个采样点，Np=pj|pj∈NP是pi的k邻近点集(本文k取经验值为15)。Np的3×3协方差矩阵C3×3定义如下

式中，为Np点集的重心；|Np|是Np点集中点的数量。

通过特征分解计算协方差矩阵C3×3的特征值λ0、λ1、λ2(0≤λ0≤λ1≤λ2)，点pi的平整度fp定义如下

fp越小，表明pi更可能位于平面上，也就是更可能属于建筑物点。计算完所有非地面点的特征fp后，每一个DSM格网的平整度取格网内所有LiDAR点的fp均值，对于没有对应LiDAR点的格网采用1.1节所提的内插方法进行内插。为保证效率，只对Mng中值为“1”的格网进行内插。

点云平整度fp描述了局部平面特征，但是当采样点位于建筑物屋脊处或者屋顶面比较复杂时，fp将会更倾向于指示为植被区域。因此本文采用另一个局部几何特征——法向量分布。通常情况下，建筑物区域局部点云的法向量会集中于某一个或某几个方向，而植被区域局部点云的法向量更发散。本文提出基于法向量方向分布直方图统计的平方变异系数来描述这一特征。令Npn是pi的kn邻近点集，首先计算Npn中所有点的法向量与竖直方向的夹角αV(αV∈[0, 90°])，然后根据αV将所有点统计到不同的直方图区间。如果pi位于植被区域，则Npn中的点更倾向于均匀分布到不同的直方图区间，而建筑物反之。因此采用此直方图的平方变异系数描述法向量分布特征fN

式中

式中，ndi是直方图区间数量；ni(i=1, 2, …, ndi)是每一个区间的点的数量；|Npn|是点集Npn中点的数量。fN越小，表明Npn中的点更倾向于均匀分布，即更可能属于植被区域，反之，则为建筑物区域。对于fN，需要设置的参数为邻域点数量kn以及直方图区间数量ndi，本文将kn设置为60，ndi设置为6。计算完所有非地面点的特征fN后，与平整度一样，内插每一个DSM格网的fN。

fp和fN是基于点的特征，通过点的邻域计算，所以Np和Npn中的点可能来自于不同的地物，比如建筑物和与其靠近的植被。此外，当点云密度较低时，对于面积小的建筑物，采样邻域点时很有可能会采样到其他地物的LiDAR点。因此，本文引入基于格网的特征来补充基于点的特征。植被区域高程分布复杂，因此nDSM生成的影像(简称归一化高程影像)中表现为丰富的纹理，而建筑物区域高程分布单一，高程影像表现纹理单一，因此本文引入灰度共生矩阵一致性来描述这一纹理特征fTH(见文献[23])。fTH越小，表明一致性越差，也就是纹理越丰富，指示该格网更可能属于植被区域。为了保证不同区域的阈值一致，本文将待提取区域的高差归一化到0~ 60 m，也就是归一化高程高于60 m的地物直接认为是建筑物，这样对于所有区域，每一个灰度级代表的高程值一致。为了加快计算速度，只对Mng中值为“1”的格网计算fTH。

影像具有丰富的颜色信息，对于建筑物与植被具有较好的区分度。本文采用NDVI来区分植被与建筑物，定义如下

式中，IR、R分别为对应颜色波段的值，采用阴影增强后的影像计算。计算完正射影像的NDVI后，将其赋予每个DSM格网。此外，生成一张基于NDVI的植被区域掩膜Mv，将fNDVI>0.3的格网值设为“1”，表示植被区域，其他区域设为“0”表示非植被区域。

以上介绍的4种特征fp、fN、fTH、fNDVI取值范围不同，因此利用Logistic方程进行归一化，定义如下

式中，x0表示归一化的中心；k控制Logistic方程的深度；k取值对结果影响较小，根据方程曲线的形态将4种特征分别设置为-35、2.0、0.2、-10；x0通过试验进行确定，即单独利用某一个特征进行建筑物提取，通过调整x0值，提取结果最优的为该特征的x0取值，对于4个特征，分别设置为0.06、0.8、18、0.12。

上节介绍的4个特征仅仅描述一个点或者一个格网的特征，并没有考虑上下文关系，因此本文引入图割算法来考虑邻域关系以保证邻域内的一致性。图割算法基本思想是构建一个权重图，每一条边的权表示相应的能量值，然后利用最大流/最小割算法寻找图的最小割[17]。图割算法的目的是通过最小化如式(7)定义的能量函数为每一个结点赋予一个标记lp。

式中，为数据项；为平滑项；β为平滑项系数；Dp(lp)表示标记lp符合结点p的程度；Vp, q(lp, lq)表示结点p和q的相似程度。

如式(7)所示，能量函数由数据项和平滑项构成，本文建筑物提取的结果是将DSM格网分为两类，一类是建筑物(building)，其他为非建筑物(non-building)，即lp∈{building, non-building}，因此数据项定义如下：

式中，λ1、λ2分别为几何和颜色特征的权重，λ11、λ12、λ13分别为几何特征中fp、fN、fTH的权重，且满足λ1+λ2=1、λ11+λ12+λ13=1。对于Mng中值为“0”的格网，直接赋予标记为“non-building”，值为“1”的格网，根据式(8)设置数据项。

由于DSM格网高程值在局部区域内具有较高的一致性，而在地物边缘处则有较明显的差异，尤其是建筑物边缘，因此本文采用格网高程值来构建平滑项。此外，在建筑物附近有高程相近的植被的情况下，可能会出现过度传播现象，易出现植被误检为建筑物或者建筑物被漏检的问题。因此本文在平滑项中引入基于NDVI指数的惩罚因子δ，平滑项设置如下

式中，hp和hq分别为相邻格网的高程值；常数ε保证分母不为0，可由LiDAR数据的高程精度确定，本文设为0.2 m；δ由植被区域掩膜Mv确定是否给予平滑项，当相邻格网在Mv中的值不同时，给予惩罚，避免分为一类，本文δ取值0.1；平滑项系数β与建筑物提取数据(如建筑物高度、建筑物屋顶复杂度等)有关，本文将其设置为1。

能量函数定义后，采用标准最大流/最小割算法[24]求解能量函数的最小割，然后确定每一个DSM格网的标记lp。

上一节图割算法优化后的结果中可能存在一些低矮地物，如灌木、篱笆等，因此进一步采用高程阈值Th和面积阈值Ta去除这些地物，即高程小于Th的格网直接去除，然后采用连通分析将提取结果聚为不同的对象，当对象面积小于Ta时直接去除。本文Th和Ta分别设置为1.5 m和2.5 m2。

经过图割优化后提取的建筑物结果可以视为一张二值图。在影像信息的辅助下，本文引入前后景分割的思想，利用图割算法对建筑物的边缘提取结果进一步优化，能量函数依然采用式(7)，但数据项和平滑项的设置与初始建筑物提取阶段有所区别。首先，将初始检测的建筑物区域设为前景，非建筑物区域设为后景，然后检测初始建筑物区域的边缘像素，并在边缘像素附近划定缓冲区。通过图割优化，将前景和后景的标记传播到缓冲区，从而实现建筑物边缘的精确提取。3个区域的数据项分别设置如下

平滑项由对应正射影像颜色信息构建，定义如下

式中，R、G、IR分别为对应颜色波段的值；εcol设置为0.1。定义能量函数后，采用标准最大流/最小割算法[24]解算能量函数对缓冲区域的格网重新标记。从能量函数的数据项可以看出，本文将初始提取的建筑物区域以及非建筑物区域设置为固定的值，防止在能量最小化过程中标记发生变化，而主要依赖于平滑项将建筑物区域和非建筑区域的标记传播到缓冲区，从而实现缓冲区的精确标记。在此基础上，进一步利用形态学滤波来消除建筑物边缘的锯齿现象。

为了验证本文方法，采用如图 2所示的ISPRS Vaihingen的数据进行测试，数据采集自德国的Vaihingen地区，具有正射影像和激光点云数据，正射影像如图 2(a)所示，分辨率为0.09 m，点云平均密度为4 pts/m2，内插的DSM如图 2(b)所示，覆盖面积约2.0 km2，拥有超过1800栋建筑物。首先采用如图 2(c)、(d)、(e)所示的3个区域基准数据测试，正射影像黄色线内为测试区域，参考数据由ISPRS官方提供，其中区域1为形状复杂的密集历史建筑物区域(图 2(c))，区域2为树木环绕的高耸居民楼区域(图 2(d))，区域3为带有小附属建筑物的居民区(图 2(e))。整个数据区域没有官方的参考数据，笔者采用手工标记的形式在正射影像上标记了所有的房屋区域，以用于验证本文方法的有效性。

本文中的λ1、λ2以及λ11、λ12、λ13通过试验确定，分别取λ1=λ2=0.5、λ11=0.25、λ12=0.5、λ13=0.25。对于ISPRS Vaihingen基准数据和全部数据，采用相同的参数设置。

图 3列举了区域1边缘优化前后建筑物区域在正射影像上套合的结果。图中右上角的区域为图中心黄色方框内区域放大显示结果。从图中可以看出，经过优化后，边缘更准确，消除了锯齿现象，并在一定程度上恢复了拐角信息。

为了进一步评价本文方法，采用ISPRS官方提供的参考数据进行定量评价。采用的指标包括基于面积(per-area)、目标(per-object)、面积大于50 m2的目标(per-object>50 m2)的完整率(%)，正确率(%)以及质量(%)(见文献[25])。结果如表 1所示，与ISPRS测试网站上其他方法结果的对比如表 2所示。

从表 1和表 2可以看出，本文方法具有一定的优势。对于基于面积的指标，本文方法取得了最优的完整率，质量指标最高。基于目标的指标中，正确率为100%，说明本文方法没有误检目标，且完整率和质量指标也较优异。基于目标面积大于50 m2的指标，全部达到了100%。此外，通过表 1可以看出，本文方法对于3个环境差异较大的区域都取得了较好的提取结果，也说明了本文方法的稳定性。基于面积的评价结果如图 4所示，黄色为正确提取，红色为误检，蓝色为漏检。从视觉上看，本文方法取得了较好的结果。误检主要来自于建筑物附近的植被区域，主要是因为这些区域处于阴影之中，即使进行了阴影增强，也难以有效区分。漏检建筑物主要为较低矮的建筑物，漏检主要原因为：有些建筑物较矮小，只含有少量激光点而无法被提取，有些建筑物部分被植物遮挡而无法提取。此外，一些屋顶结构复杂的建筑物以及在地面点滤波时错分为地面的建筑物也被漏检。本文优化建筑物边缘信息时，采用的是正射影像，因此，正射影像建筑物边缘纠正的效果也对本文方法的结果产生一定影响。

由于ISPRS Vaihingen全区域没有提供参考数据，因此笔者通过人机交互方式进行了标记。对于Vaihingen全区域测试了基于面积(per-area)、基于目标(per-object)两个指标，结果如表 3所示，基于面积的评价结果如图 5所示。

从表 3和图 5可以看出，本文方法在全部大范围数据上也表现出了较好的性能，说明本文方法具有一定的适应性。从图 5中还可以看出，主要的漏检建筑物面积较小，另外部分建筑物屋顶覆盖了较多的植被也被漏检。除了部分误检是由植被区域导致以外，还有部分误检是由于地面点在点云滤波时被错分为非地面，这部分区域被误检为建筑物。该区域右下角为工厂，有大量货物堆积在空旷地面上，这些货物也被误检为建筑物。

本文方法采用C++实现，运行环境为台式机(Intel i5-3470CPU，主频：3.20 GHz，内存8 GB)，对于3个测试区域的提取时间分别为12、15和14 s，全部区域的提取时间为660 s，其中耗时最多的步骤为点云几何特征计算。总体上看，本文方法具有较好的效率。

本文基于图割算法，提出融合LiDAR点云与正射影像的建筑物自动提取方法。本文方法充分利用LiDAR点云与影像提供的几何以及颜色特征，通过图割算法进行优化，获得了邻域平滑一致的结果。综合DSM和NDVI构建能量函数平滑项，有效防止了建筑物与植被区域的标记过度传播。利用影像分辨率高的特点，基于前后景分割的方法进一步优化建筑物边缘，获得了更加准确的边缘信息。通过ISPRS Vaihingen地区的数据进行测试，并与其他方法进行比较，证实了本文方法可以有效地提取建筑物，精度较高，具有一定的实用推广价值。进一步顾及建筑物边缘的规则特性提高建筑物的检测精度是下一步的研究内容。